use "[path to file here].dta"



***/// setup 
xtset id_cow year
sort id_cow year

***/// generating study variables

	*coutry-year counter
by id_cow: generate time = _n

	*rescaling to range between 0 and 100
gen rev_source = (state_rev_source+3.379)*10

	*rescaling democracy index to range between 0 and 100
gen elect_recode = elect_dem_vdem*100

	*creating IGF dummies
gen igf_recode= int_forum_igf
replace igf_recode =2 if year<2006 & int_forum_igf !=.

	*creating fdi_inflow pct of gdp measure
gen fdi_raw= fdi_inward_mil_curr_dol*1000000
gen fdi_pct_gdp = (fdi_raw/gdp)*100

	*creating rescaled export concentration index
gen expt_concentration_recode = expt_concentration_index*100
sum expt_concentration_recode expt_concentration_index
drop expt_concentration_index
sum expt_

	*creating debt as a % GDP measure (for robusteness checks below)
gen debt_pct_gdp = (tot_ext_debt_curr_us/gdp)*100





*****/////  Model to designate sample (model 6)  /////*****

*get alpha value
nbreg mob_tel i.inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 merch_ ///
rev_source rur_pop_pct ngo fdi_pct_ expt_c ib2.igf_r ///
c.time##c.time if inc_grp95!=. & inc_grp95!=4, ///
vce(robust) exposure(pop_tot)



*dev sample
xtgee mob_tel i.inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 merch_ ///
rev_source rur_pop_pct ngo fdi_pct_ expt_c ib2.igf_r ///
c.time##c.time if inc_grp95!=. & inc_grp95!=4, ///
family(nbinomial 1.168) link(log) ///
corr(ar1) vce(robust) exposure(pop_tot) eform force

gen dev_samp = 0
replace dev_samp = 1 if e(sample)==1





****************************** SUMMARY STATS **********************************

	*creating dummies (g1-g4 and IGF) and quadratic for proceedures that don't take operators
tabulate inc_grp95, gen(g)
tabulate igf_, gen(igf)
gen time2 = time*time



* Table 1 - betwen and within panel summary stats
gen mil_mob=mob_tel/1000000

xtsum mil_mob gdp_grow elect_r work_pop_chg fixed_tel_100 merch_ ///
rev_source ngo rur_pop_pct fdi_pct_ expt_c if dev_sam==1



* Appendix 1 - spearman correlations 
spearman mob_tel gdp_grow elect_r work_pop_chg fixed_tel_100 merch ///
rev_source ngo rur_pop_pct fdi_pct_ expt_c time pop_tot if dev_sam==1





**************************** residual analyses  ********************************

******************** assesing residuals BETWEEN countries **********************

****/// Pearson (standardized) observation residuals and predicted values
****/// and averaged by panel - need to remove Djibouti here 
****/// (excluded with AR1 in final models)

*model 6 (ind corr and no robust to test naive DGP)
xtgee mob_tel i.inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 merch_ ///
rev_source rur_pop_pct ngo fdi_pct_ expt_c ib2.igf_r ///
c.time##c.time if inc_grp95!=. & inc_grp95!=4 & id_!=522, ///
family(nbinomial 1.168) link(log) ///
corr(ind) exposure(pop_tot) eform force



*Pearson residuals by observation
predict yhatm if e(sample)==1, rate
gen residm = mob_tel-yhatm if e(sample)==1
egen sd_residm=sd(residm) if e(sample)==1
gen std_resid_o = residm/sd_residm if e(sample)==1



*Standardized predicted value by observation
sum yhatm if e(sample)==1, meanonly
gen yhatm_diff = yhatm - r(mean) if e(sample)==1
egen sd_yhatm=sd(yhatm) if e(sample)==1
gen std_yhatm_o = yhatm_diff/sd_yhatm if e(sample)==1

twoway (scatter std_resid_o yhatm if year==1995, mlabel(iso)) ///
(lfit std_resid_o yhatm if year==1995)
graph save 1995, replace

twoway (scatter std_resid_o yhatm if year==2000, mlabel(iso)) ///
(lfit std_resid_o yhatm if year==2000)
graph save 2000, replace

twoway (scatter std_resid_o yhatm if year==2005, mlabel(iso)) ///
(lfit std_resid_o yhatm if year==2005)
graph save 2005, replace

twoway (scatter std_resid_o yhatm if year==2010, mlabel(iso)) ///
(lfit std_resid_o yhatm if year==2010)
graph save 2010, replace

twoway (scatter std_resid_o yhatm if year==2014, mlabel(iso)) ///
(lfit std_resid_o yhatm if year==2014)
graph save 2014, replace

graph combine 1995.gph 2000.gph 2005.gph 2010.gph 2014.gph



*average pearson residuals by panel
sort id_ year
by id_: egen avg_std_resid_o=mean(std_resid_o) if e(sample)==1

*same for standardized pred values
sort id_ year
by id_: egen avg_std_yhatm_o=mean(std_yhatm_o) if e(sample)==1



*inspect mean pearson resid by panel (larger than absolute value of 2 is cutoff),
scatter avg_std_resid_o id_c if e(sample)==1, mlabel(iso)

*mean pearson vs. mean fitted
twoway (scatter avg_std_resid_o avg_std_yhatm_o if e(sample)==1, mlabel(iso)) ///
(lfit avg_std_resid_o avg_std_yhatm_o if e(sample)==1, mlabel(iso))

*URY, CHN, IND identified as potentially problematic







********************* assesing residuals WITHIN countries **********************


*model 5 (using ind corr and no robust to test naive DGP)
xtgee mob_tel i.inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 merch_ ///
rev_source rur_pop_pct ngo fdi_pct_ expt_c ib2.igf_r ///
c.time##c.time if inc_grp95!=. & inc_grp95!=4 & id_!=522, ///
family(nbinomial 1.168) link(log) ///
corr(ind) exposure(pop_tot) eform force

*using rate specification residuals
predict yhatr if e(sample)==1, rate
gen residr = mob_tel - yhatr



*inspection of country specific residuals over time shows clear serial 
*correlation
twoway (scatter residr time if id<200 & id!=70 &id!=140 &id!=160 & ///
residr!=., by(iso3))

twoway (scatter residr time if id<400 & id>=200 & id!=365 & ///
residr!=., by(iso3))

twoway (scatter residr time if id<600 & id>=400 &id!=475 & ///
residr!=., by(iso3))

twoway (scatter residr time if id<951 & id>=600 & id!=710 & id!=750 & ///
id!=850 & residr!=., by(iso3))




****//// Runs Test - Non-parametric test of random errors (raw residuals)

*see Hardin and Hilbe (2013) pg 172 and Chang (2000)

*np=# positive; nn=# neg, T=# of runs
*E(T) = ((2*np*nn)/(np+nn))+1
*V(T) = ((2*np*nn)*((2*np*nn)-np-nn))/(((np+nn)^2)*(np+nn-1))
*Wz(t-stat) = (T-E(T))/sqrt(V(T))


gen sign = .
replace sign = 1 if (residr >0) & (residr!=.)
replace sign = 0 if (residr ==0)
replace sign = -1 if (residr <0) & (residr!=.)

*identify number of positive and negative resid values (nn and np)
tab sign
sort id time

*identify number of spells
tsspell sign
egen nspells = max(_spell) if e(sample)==1, by(id)
tab nspells

*copy/paste table in excel and sum ONLY the number of runs attached to an id 
*(ignore the "Total" value at the bottom of the table) to get T (# of runs)
table id_cow, statistic(max nspells)

****/// formula from Hardin and Hilbe (2013:172) and Chang (2000)

*np=2379; nn=95, T=210

*E(T) = ((2*np*nn)/(np+nn))+1
di ((2*2379*95)/(2379+95))+1
*=183.7

*V(T) = ((2*np*nn)*((2*np*nn)-np-nn))/(((np+nn)^2)*(np+nn-1))
di ((2*2379*95)*((2*2379*95)-2379-95))/(((2379+95)^2)*(2379+95-1))
*=13.42

*Wz(t-stat) = (T-E(T))/sqrt(V(T))
di (210-183.7)/sqrt(13.42)
*= 7.18 - reject null of independence (randomness) in residuals, 
*in conjunction with scatter plot of within country residuals, the runs test 
*indicates some form of dependence (probably autocorelation)




*****///// Testing specifically for serial correlation

tsset id time


*Drukker (2003)
xtserial mob_tel gdp_grow elect_r fixed_tel_100 merch_ ///
rev_source rur_pop_pct ngo pop_tot work_pop_chg ///
fdi_pct_ expt_c igf1 igf2 if dev_samp==1, output 
	*conclude serial corr
	
xtserial mob_tel gdp_grow elect_r fixed_tel_100 merch_ ///
rev_source rur_pop_pct ngo pop_tot work_pop_chg time time2 ///
fdi_pct_ expt_c igf1 igf2 if dev_samp==1, output
	*conclude serial corr

	



********************************* influence ***********************************


*jackknife DFBETA (Canette 2014, Stata blog)

local vars g2 g3 gdp_grow elect_r fixed_tel_100 merch_ ///
rev_source rur_pop_pct ngo work_pop_chg fdi_pct_ expt_c igf1 igf2 time time2

gen y = mob_tel       // short name

jackknife _b _se, cluster(id_cow) keep: xtgee y `vars' if dev_samp==1, ///
family(nbinomial 1.168) link(log) corr(ar1) exposure(pop_tot) vce(robust) force

sum _*  //pseudovalues

xtgee y `vars' if dev_samp==1, family(nbinomial 1.168) link(log) ///
corr(ar1) exposure(pop_tot) vce(robust) force //entire sample

foreach z in `vars'{
gen double b_`z' = (e(N)*_b[`z'] - _b_`z')/(e(N) -1)
gen double se_`z' = (e(N)*_se[`z'] - _se_`z')/(e(N)-1)
gen double dfb_`z' = (_b[`z']-b_`z')/se_`z'
label var dfb_`z' "DFBETAS `z'"
}


* cuttoff is +/- 0.04 (2/sqrt of n or 2/sqrt(2474)) but using conservative ~0.03
scatter dfb_g2 id_, mlabel(iso3) // THA high
scatter dfb_g3 id_, mlabel(iso3) // 
scatter dfb_gdp id_, mlabel(iso3) // THA high, LBY low
scatter dfb_elect id_, mlabel(iso3) // 
scatter dfb_fix id_, mlabel(iso3) // 
scatter dfb_merch_ id_, mlabel(iso3) // MYS high
scatter dfb_rev id_, mlabel(iso3) //THA low
scatter dfb_rur id_, mlabel(iso3) // THA high
scatter dfb_ngo id_, mlabel(iso3) // THA high
scatter dfb_work id_, mlabel(iso3) //
scatter dfb_fdi_pct_ id_, mlabel(iso3) // MLT low
scatter dfb_expt id_, mlabel(iso3) //
scatter dfb_igf1 id_, mlabel(iso3) // 
scatter dfb_igf2 id_, mlabel(iso3) //
scatter dfb_time id_, mlabel(iso3) // THA low
scatter dfb_time2 id_, mlabel(iso3) // THA high, of course

*THA, LBY, MYS, and MLT potentially influential and
*CHN, IND, URY from residual analysis above

list id_ if iso3=="THA" // 800
list id_ if iso3=="LBY" // 620
list id_ if iso3=="MYS" // 820
list id_ if iso3=="URY" // 165
list id_ if iso3=="CHN" // 710
list id_ if iso3=="IND" // 750
list id_ if iso3=="MLT" // 338





*re-estimate for sensitivity analyses

*model 6 from Table 2 (for comparison)
xtgee mob_tel i.inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 merch_ ///
rur_pop_pct rev_source ngo fdi_pct_ expt_c ib2.igf_r ///
c.time##c.time if dev_sam==1, family(nbinomial 1.168) link(log) ///
corr(ar1) vce(robust) exposure(pop_tot) eform force

*model 6 without Libya - elect to marg sig
xtgee mob_tel i.inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 merch_ ///
rur_pop_pct rev_source ngo fdi_pct_ expt_c ib2.igf_r ///
c.time##c.time if dev_sam==1 & id_!=620, family(nbinomial 1.168) link(log) ///
corr(ar1) vce(robust) exposure(pop_tot) eform force

*model 6 without Malaysia - no change
xtgee mob_tel i.inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 merch_ ///
rur_pop_pct rev_source ngo fdi_pct_ expt_c ib2.igf_r ///
c.time##c.time if dev_sam==1 & id_!=820, family(nbinomial 1.168) link(log) ///
corr(ar1) vce(robust) exposure(pop_tot) eform force

*model 6 without Thailand - no change
xtgee mob_tel i.inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 merch_ ///
rur_pop_pct rev_source ngo fdi_pct_ expt_c ib2.igf_r ///
c.time##c.time if dev_sam==1 & id_!=800, family(nbinomial 1.168) link(log) ///
corr(ar1) vce(robust) exposure(pop_tot) eform force

*model 6 without Uruguay - no change
xtgee mob_tel i.inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 merch_ ///
rur_pop_pct rev_source ngo fdi_pct_ expt_c ib2.igf_r ///
c.time##c.time if dev_sam==1 & id_!=165, family(nbinomial 1.168) link(log) ///
corr(ar1) vce(robust) exposure(pop_tot) eform force

*model 6 without India - no change
xtgee mob_tel i.inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 merch_ ///
rur_pop_pct rev_source ngo fdi_pct_ expt_c ib2.igf_r ///
c.time##c.time if dev_sam==1 & id_!=750, family(nbinomial 1.168) link(log) ///
corr(ar1) vce(robust) exposure(pop_tot) eform force

*model 6 without China - no change
xtgee mob_tel i.inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 merch_ ///
rur_pop_pct rev_source ngo fdi_pct_ expt_c ib2.igf_r ///
c.time##c.time if dev_sam==1 & id_!=710, family(nbinomial 1.168) link(log) ///
corr(ar1) vce(robust) exposure(pop_tot) eform force

*model 6 without Malta - fdi drops from sig
xtgee mob_tel i.inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 merch_ ///
rur_pop_pct rev_source ngo fdi_pct_ expt_c ib2.igf_r ///
c.time##c.time if dev_sam==1 & id_!=338, family(nbinomial 1.168) link(log) ///
corr(ar1) vce(robust) exposure(pop_tot) eform force


*IN NEARLY ALL CASES, NO CHANGE (Malta is an odd case in FDI in reality)

*model 6 all excluded - Some movement here, but mostly shifting between conventional and marginal significance. And the focal variables (rurality, revenue source) and the majority of all variables remain consistent. 
xtgee mob_tel i.inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 merch_ ///
rur_pop_pct rev_source ngo fdi_pct_ expt_c ib2.igf_r ///
c.time##c.time if dev_sam==1 & id_!=620 & id_!=820 ///
& id_!=800 & id_!=165 & id_!=750 & id!=710 & id!=338, family(nbinomial 1.168) ///
link(log) corr(ar1) vce(robust) exposure(pop_tot) eform force





********************* specifying the correlation structure *********************

*model 6 (Table 2) (ind corr)
xtgee mob_tel i.inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 merch_ ///
rur_pop_pct rev_source ngo fdi_pct_ expt_c ib2.igf_r ///
c.time##c.time if dev_sam==1, family(nbinomial 1.168) link(log) ///
corr(ind) vce(robust) exposure(pop_tot) eform force

*blank matrix, as assumption is independent
estat wcorr

*qic value
qic mob_tel g2 g3 gdp_grow elect_r work_pop_chg fixed_tel_100 merch_ ///
rur_pop_pct rev_source ngo fdi_pct_ expt_c igf1 igf2 ///
time time2 if dev_sam==1, family(nbinomial 1.168) link(log) ///
corr(ind) robust exposure(pop_tot) eform force i(id_cow) t(year) nodis

*=5743



*model 6 (Table 2) (AR1)
xtgee mob_tel i.inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 merch_ ///
rur_pop_pct rev_source ngo fdi_pct_ expt_c ib2.igf_r ///
c.time##c.time if dev_sam==1, family(nbinomial 1.168) link(log) ///
corr(ar1) vce(robust) exposure(pop_tot) eform force

estat wcorr
*clear pattern of AR correlation in residuals (decaying over time)

*qic value
qic mob_tel g2 g3 gdp_grow elect_r work_pop_chg fixed_tel_100 merch_ ///
rur_pop_pct rev_source ngo fdi_pct_ expt_c igf1 igf2 ///
time time2 if dev_sam==1, family(nbinomial 1.168) link(log) ///
corr(ar1) robust exposure(pop_tot) eform force i(id_cow) t(year) nodis

*=7123



















****************/////////////// MODELS (TABLE 2) \\\\\\\\\\\\\\\**************** 



**************************** model 1

*get alpha value
nbreg mob_tel i.inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 ///
rur_pop_pct c.time##c.time if dev_sam==1, ///
vce(robust) exposure(pop_tot) 



*model 1
xtgee mob_tel i.inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 ///
rur_pop_pct c.time##c.time if dev_sam==1, family(nb 1.224) link(log) ///
corr(ar1) vce(robust) exposure(pop_tot) eform force

testnl (_b[2.inc_grp95]=0) (_b[3.inc_grp95]=0)



*qic value
qic mob_tel g2 g3 gdp_grow elect_r work_pop_chg fixed_tel_100 ///
rur_pop_pct time time2 if dev_sam==1, family(nbinomial 1.224) link(log) ///
corr(ar1) robust exposure(pop_tot) eform force i(id_cow) t(year) nodis



*collinearity - inc_grp dummies and pre_igf dummy excluded
collin rur_pop_pct inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 time if dev_sam==1





************************** model 2

*get alpha value
nbreg mob_tel i.inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 rur_pop_pct merch_ ///
rev_source fdi_pct_ ///
c.time##c.time if dev_sam==1, ///
vce(robust) exposure(pop_tot)



*model 2
xtgee mob_tel i.inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 rur_pop_pct merch_ ///
 rev_source fdi_pct_ ///
c.time##c.time if dev_sam==1, family(nbinomial 1.181) link(log) ///
corr(ar1) vce(robust) exposure(pop_tot) eform force

testnl (_b[2.inc_grp95]=0) (_b[3.inc_grp95]=0)



*qic value
qic mob_tel g2 g3 gdp_grow elect_r work_pop_chg fixed_tel_100 rur_pop_pct merch_ ///
rev_source fdi_pct_ ///
time time2 if dev_sam==1, family(nbinomial 1.181) link(log) ///
corr(ar1) robust exposure(pop_tot) eform force i(id_cow) t(year) nodis



*collinearity - inc_grp dummies and pre_igf dummy excluded
collin rur_pop_pct inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 ///
rev_source merch_ fdi_pct_ ///
time if dev_sam==1





************************** model 3

*get alpha value
nbreg mob_tel i.inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 rur_pop_pct ///
rev_source expt_c ///
c.time##c.time if dev_sam==1, ///
vce(robust) exposure(pop_tot)



*model 3
xtgee mob_tel i.inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 rur_pop_pct ///
rev_source fdi_pct_ expt_c ///
c.time##c.time if dev_sam==1, family(nbinomial 1.196) link(log) ///
corr(ar1) vce(robust) exposure(pop_tot) eform force

testnl (_b[2.inc_grp95]=0) (_b[3.inc_grp95]=0)



*qic value
qic mob_tel g2 g3 gdp_grow elect_r work_pop_chg fixed_tel_100 rur_pop_pct ///
rev_source fdi_pct_ expt_c ///
time time2 if dev_sam==1, family(nbinomial 1.196) link(log) ///
corr(ar1) robust exposure(pop_tot) eform force i(id_cow) t(year) nodis



*collinearity - inc_grp dummies and pre_igf dummy excluded
collin rur_pop_pct inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 ///
rev_source fdi_pct_ expt_c ///
time if dev_sam==1





************************** model 4

*get alpha value
nbreg mob_tel i.inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 rur_pop_pct merch_ ///
rev_source fdi_pct_ expt_c ///
c.time##c.time if dev_sam==1, ///
vce(robust) exposure(pop_tot)



*model 4
xtgee mob_tel i.inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 rur_pop_pct merch_ ///
rev_source fdi_pct_ expt_c ///
c.time##c.time if dev_sam==1, family(nbinomial 1.17) link(log) ///
corr(ar1) vce(robust) exposure(pop_tot) eform force

testnl (_b[2.inc_grp95]=0) (_b[3.inc_grp95]=0)



*qic value
qic mob_tel g2 g3 gdp_grow elect_r work_pop_chg fixed_tel_100 ///
rur_pop_pct merch_ ///
rev_source fdi_pct_ expt_c ///
time time2 if dev_sam==1, family(nbinomial 1.17) link(log) ///
corr(ar1) robust exposure(pop_tot) eform force i(id_cow) t(year) nodis



*collinearity - inc_grp dummies and pre_igf dummy excluded
collin rur_pop_pct inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 ///
rev_source merch_ fdi_pct_ expt_c ///
time if dev_sam==1





************************** model 5

*get alpha value
nbreg mob_tel i.inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 ngo ///
ib2.igf_r rur_pop_pct ///
c.time##c.time if dev_sam==1, ///
vce(robust) exposure(pop_tot)



*model 5
xtgee mob_tel i.inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100  rur_pop_pct ///
ngo ib2.igf_r  ///
c.time##c.time if dev_sam==1, family(nbinomial 1.223) link(log) ///
corr(ar1) vce(robust) exposure(pop_tot) eform force

testnl (_b[2.inc_grp95]=0) (_b[3.inc_grp95]=0)
testnl (_b[0.igf_recode]=0) (_b[1.igf_recode]=0)



*qic value
qic mob_tel g2 g3 gdp_grow elect_r work_pop_chg fixed_tel_100 rur_pop_pct ///
ngo igf1 igf2 ///
time time2 if dev_sam==1, family(nbinomial 1.223) link(log) ///
corr(ar1) robust exposure(pop_tot) eform force i(id_cow) t(year) nodis



*collinearity - inc_grp dummies and pre_igf dummy excluded
collin inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 ///
rur_pop_pct ngo igf_r ///
time if dev_sam==1





************************** model 6

*get alpha value
nbreg mob_tel i.inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 rur_pop_pct merch_ ///
rev_source fdi_pct_ expt_c ngo ib2.igf_r ///
c.time##c.time if dev_sam==1, ///
vce(robust) exposure(pop_tot)



*model 6
xtgee mob_tel i.inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 rur_pop_pct merch_ ///
rev_source fdi_pct_ expt_c ngo ib2.igf_r ///
c.time##c.time if dev_sam==1, family(nbinomial 1.168) link(log) ///
corr(ar1) vce(robust) exposure(pop_tot) eform force

testnl (_b[2.inc_grp95]=0) (_b[3.inc_grp95]=0)
testnl (_b[0.igf_recode]=0) (_b[1.igf_recode]=0)



*qic value
qic mob_tel g2 g3 gdp_grow elect_r work_pop_chg fixed_tel_100 rur_pop_pct merch_ ///
rev_source fdi_pct_ expt_c ngo igf1 igf2 ///
time time2 if dev_sam==1, family(nbinomial 1.168) link(log) ///
corr(ar1) robust exposure(pop_tot) eform force i(id_cow) t(year) nodis



*collinearity - inc_grp dummies and pre_igf dummy excluded
collin inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 ///
rur_pop_pct rev_source ngo merch_ fdi_pct_ expt_c igf_r ///
time if dev_sam==1





*********************** RURALITYxREV SOURCE interaction ***********************




***************************** marginal effects *********************************

*Model 7 - pop_ruralXrev_source interaction
xtgee mob_tel i.inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 merch_ ///
fdi_pct_ expt_c ngo ib2.igf_r c.rev_source##c.rur_pop_pct ///
c.time##c.time if dev_sam==1, family(nbinomial 1.168) link(log) ///
corr(ar1) vce(robust) exposure(pop_tot) force eform

sum rev_source if dev_s==1, d
sum rur_pop_pct if dev==1, d



*at 10th and 90th pct
margins, at(rev_source =(25 40 55) rur_pop_pct=(0(10)100) (means) _all) vsquish
mplotoffset, x(rur_pop_pct) offset(1.1) recastci(rspike)
*graph save new_rev, replace



*qic value
gen revXrur = rev_source*rur_pop_pct

qic mob_tel g2 g3 gdp_grow elect_r work_pop_chg fixed_tel_100 merch_ ///
rur_pop_pct rev_source fdi_pct_ expt_c ngo igf1 igf2 revXrur ///
time time2 if dev_sam==1, family(nbinomial 1.168) link(log) ///
corr(ar1) robust exposure(pop_tot) eform force i(id_cow) t(year) nodis 




********************** Test using multiplicative effects ***********************

*Using procedure described by Hilbe (2011) pp.528-529

*constructing rural decile dummies
egen rur_dec=xtile(rur_pop_pct) if dev_s==1, nq(10)
tab rur_dec if dev_s==1, gen(rur) 

*costructing decileXrev_source interactions (1st decile is ref)
gen revXrur2 = rur2*rev_source
gen revXrur3 = rur3*rev_source
gen revXrur4 = rur4*rev_source
gen revXrur5 = rur5*rev_source
gen revXrur6 = rur6*rev_source
gen revXrur7 = rur7*rev_source
gen revXrur8 = rur8*rev_source
gen revXrur9 = rur9*rev_source
gen revXrur10 = rur10*rev_source

*running model 7 with decile interactions
xtgee mob_tel i.inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 merch_ ///
fdi_pct_ expt_c ngo ib2.igf_r ///
rev_source rur2 rur3 rur4 rur5 rur6 rur7 rur8 rur9 rur10 revXrur2 revXrur3 revXrur4 revXrur5 revXrur6 revXrur7 revXrur8 revXrur9 revXrur10 ///
c.time##c.time if dev_sam==1, family(nbinomial 1.168) link(log) ///
corr(ar1) vce(robust) exposure(pop_tot) eform force

*list model vectors
ereturn list
matrix list e(b)

sum rev_source if dev_s==1, d
***/// calculate multiplicative values for interactions 
***/// IRR=exp[b_dum+(b_int*rev_source_value)]

***/// calculate CIs
***/// SE= sqrt(V(b_dum)+((rev_val^2)*V(b_int))+(2*rev_val*COV(b_dum,b_int)))
mat list e(V)
matrix a = e(V)
svmat a


**********////////// For figure 3 \\\\\\\\\\*********



********** rur level 2 (2nd most urban)
*revXrur2 - 10thpctile (value of 25):
di exp(_b[rur2]+(_b[revXrur2]*25)) //= 0.527
*revXrur2 - 50thpctile (value of 40):
di exp(_b[rur2]+(_b[revXrur2]*40)) //= 0.815
*revXrur2 - 90thpctile (value of 55):
di exp(_b[rur2]+(_b[revXrur2]*55)) //= 1.260

*revXrur2 - 10thpctile (value of 25):
display a[16,16] // var of rur2
di a[25,25] // var of revXrur2
di a[16,25] //cov of rur2, revXrur2
di sqrt(a[16,16]+((25^2)*a[25,25])+(2*25*a[16,25])) // se=0.22274145
*LB CI - use non exponentiated coefficient
di exp((_b[rur2]+(_b[revXrur2]*25)) - 1.96*0.22274145) // =0.340
*UB CI - use non exponentiated coefficient
di exp((_b[rur2]+(_b[revXrur2]*25)) + 1.96*0.22274145) // =0.815

*revXrur2 - 50thpctile (value of 40):
di sqrt(a[16,16]+((40^2)*a[25,25])+(2*40*a[16,25])) // se=0.13727838
*LB CI - use non exponentiated coefficient
di exp((_b[rur2]+(_b[revXrur2]*40)) - 1.96*0.13727838) // =0.623
*UB CI - use non exponentiated coefficient
di exp((_b[rur2]+(_b[revXrur2]*40)) + 1.96*0.13727838) // =1.066

*revXrur2 - 90thpctile (value of 55):
di sqrt(a[16,16]+((55^2)*a[25,25])+(2*55*a[16,25])) // se=0.1143285
*LB CI - use non exponentiated coefficient
di exp((_b[rur2]+(_b[revXrur2]*55)) - 1.96*0.1143285) // =1.007
*UB CI - use non exponentiated coefficient
di exp((_b[rur2]+(_b[revXrur2]*55)) + 1.96*0.1143285) // =1.577





********** rur level 3 (3nd most urban)
*revXrur3 - 10thpctile (value of 25):
di exp(_b[rur3]+(_b[revXrur3]*25)) //= 0.573
*revXrur3 - 50thpctile (value of 40):
di exp(_b[rur3]+(_b[revXrur3]*40)) //= 0.858
*revXrur3 - 90thpctile (value of 55):
di exp(_b[rur3]+(_b[revXrur3]*55)) //= 1.285

*revXrur3 - 10thpctile (value of 25):
display a[17,17] // var of rur3
di a[26,26] // var of revXrur3
di a[17,26] //cov of rur3, revXrur3
di sqrt(a[17,17]+((25^2)*a[26,26])+(2*25*a[17,26])) // se=0.26971793
*LB CI - use non exponentiated coefficient
di exp((_b[rur3]+(_b[revXrur3]*25)) - 1.96*0.26971793) // =0.338
*UB CI - use non exponentiated coefficient
di exp((_b[rur3]+(_b[revXrur3]*25)) + 1.96*0.26971793) // =0.972

*revXrur3 - 50thpctile (value of 40):
di sqrt(a[17,17]+((40^2)*a[26,26])+(2*40*a[17,26])) // se=0.16998485
*LB CI - use non exponentiated coefficient
di exp((_b[rur3]+(_b[revXrur3]*40)) - 1.96*0.16998485) // =0.615
*UB CI - use non exponentiated coefficient
di exp((_b[rur3]+(_b[revXrur3]*40)) + 1.96*0.16998485) // =1.197

*revXrur3 - 90thpctile (value of 55):
di sqrt(a[17,17]+((55^2)*a[26,26])+(2*55*a[17,26])) // se=0.14455306
*LB CI - use non exponentiated coefficient
di exp((_b[rur3]+(_b[revXrur3]*55)) - 1.96*0.14455306) // =0.968
*UB CI - use non exponentiated coefficient
di exp((_b[rur3]+(_b[revXrur3]*55)) + 1.96*0.14455306) // =1.706





********** rur level 4 (4th most urban)
*revXrur4 - 10thpctile (value of 25):
di exp(_b[rur4]+(_b[revXrur4]*25)) //= 0.584
*revXrur4 - 50thpctile (value of 40):
di exp(_b[rur4]+(_b[revXrur4]*40)) //= 0.836
*revXrur4 - 90thpctile (value of 55):
di exp(_b[rur4]+(_b[revXrur4]*55)) //= 1.197

*revXrur4 - 10thpctile (value of 25):
display a[18,18] // var of rur4
di a[27,27] // var of revXrur4
di a[18,27] //cov of rur4, revXrur4
di sqrt(a[18,18]+((25^2)*a[27,27])+(2*25*a[18,27])) // se=0.32194491
*LB CI - use non exponentiated coefficient
di exp((_b[rur4]+(_b[revXrur4]*25)) - 1.96*0.32194491) // =0.311
*UB CI - use non exponentiated coefficient
di exp((_b[rur4]+(_b[revXrur4]*25)) + 1.96*0.32194491) // =1.098

*revXrur4 - 50thpctile (value of 40):
di sqrt(a[18,18]+((40^2)*a[27,27])+(2*40*a[18,27])) // se=0.19620958
*LB CI - use non exponentiated coefficient
di exp((_b[rur4]+(_b[revXrur4]*40)) - 1.96*0.19620958) // =0.569
*UB CI - use non exponentiated coefficient
di exp((_b[rur4]+(_b[revXrur4]*40)) + 1.96*0.19620958) // =1.229

*revXrur4 - 90thpctile (value of 55):
di sqrt(a[18,18]+((55^2)*a[27,27])+(2*55*a[18,27])) // se=0.17549818
*LB CI - use non exponentiated coefficient
di exp((_b[rur4]+(_b[revXrur4]*55)) - 1.96*0.17549818) // =0.848
*UB CI - use non exponentiated coefficient
di exp((_b[rur4]+(_b[revXrur4]*55)) + 1.96*0.17549818) // =1.688





********** rur level 5 (5th most urban)
*revXrur5 - 10thpctile (value of 25):
di exp(_b[rur5]+(_b[revXrur5]*25)) //= 0.593
*revXrur5 - 50thpctile (value of 40):
di exp(_b[rur5]+(_b[revXrur5]*40)) //= 0.828
*revXrur5 - 90thpctile (value of 55):
di exp(_b[rur5]+(_b[revXrur5]*55)) //= 1.158

*revXrur5 - 10thpctile (value of 25):
display a[19,19] // var of rur5
di a[28,28] // var of revXrur5
di a[19,28] //cov of rur5, revXrur5
di sqrt(a[19,19]+((25^2)*a[28,28])+(2*25*a[19,28])) // se=0.32771095
*LB CI - use non exponentiated coefficient
di exp((_b[rur5]+(_b[revXrur5]*25)) - 1.96*0.32771095) // =0.312
*UB CI - use non exponentiated coefficient
di exp((_b[rur5]+(_b[revXrur5]*25)) + 1.96*0.32771095) // =1.126

*revXrur5 - 50thpctile (value of 40):
di sqrt(a[19,19]+((40^2)*a[28,28])+(2*40*a[19,28])) // se=0.20267661
*LB CI - use non exponentiated coefficient
di exp((_b[rur5]+(_b[revXrur5]*40)) - 1.96*0.20267661) // =0.557
*UB CI - use non exponentiated coefficient
di exp((_b[rur5]+(_b[revXrur5]*40)) + 1.96*0.20267661) // =1.232

*revXrur5 - 90thpctile (value of 55):
di sqrt(a[19,19]+((55^2)*a[28,28])+(2*55*a[19,28])) // se=0.1661408
*LB CI - use non exponentiated coefficient
di exp((_b[rur5]+(_b[revXrur5]*55)) - 1.96*0.1661408) // =0.836
*UB CI - use non exponentiated coefficient
di exp((_b[rur5]+(_b[revXrur5]*55)) + 1.96*0.1661408) // =1.604





********** rur level 6 (6th most urban)
*revXrur6 - 10thpctile (value of 25):
di exp(_b[rur6]+(_b[revXrur6]*25)) //= 0.515
*revXrur6 - 50thpctile (value of 40):
di exp(_b[rur6]+(_b[revXrur6]*40)) //= 0.723
*revXrur6 - 90thpctile (value of 55):
di exp(_b[rur6]+(_b[revXrur6]*55)) //= 1.015

*revXrur6 - 10thpctile (value of 25):
display a[20,20] // var of rur6
di a[29,29] // var of revXrur6
di a[20,29] //cov of rur6, revXrur6
di sqrt(a[20,20]+((25^2)*a[29,29])+(2*25*a[20,29])) // se=0.3336413
*LB CI - use non exponentiated coefficient
di exp((_b[rur6]+(_b[revXrur6]*25)) - 1.96*0.3336413) // =0.268
*UB CI - use non exponentiated coefficient
di exp((_b[rur6]+(_b[revXrur6]*25)) + 1.96*0.3336413) // =0.991

*revXrur6 - 50thpctile (value of 40):
di sqrt(a[20,20]+((40^2)*a[29,29])+(2*40*a[20,29])) // se=0.22201576
*LB CI - use non exponentiated coefficient
di exp((_b[rur6]+(_b[revXrur6]*40)) - 1.96*0.22201576) // =0.468
*UB CI - use non exponentiated coefficient
di exp((_b[rur6]+(_b[revXrur6]*40)) + 1.96*0.22201576) // =1.117

*revXrur6 - 90thpctile (value of 55):
di sqrt(a[20,20]+((55^2)*a[29,29])+(2*55*a[20,29])) // se=0.1891319
*LB CI - use non exponentiated coefficient
di exp((_b[rur6]+(_b[revXrur6]*55)) - 1.96*0.1891319) // =0.700
*UB CI - use non exponentiated coefficient
di exp((_b[rur6]+(_b[revXrur6]*55)) + 1.96*0.1891319) // =1.470





********** rur level 7 (7th most urban)
*revXrur7 - 10thpctile (value of 25):
di exp(_b[rur7]+(_b[revXrur7]*25)) //= 0.400
*revXrur7 - 50thpctile (value of 40):
di exp(_b[rur7]+(_b[revXrur7]*40)) //= 0.671
*revXrur7 - 90thpctile (value of 55):
di exp(_b[rur7]+(_b[revXrur7]*55)) //= 1.125

*revXrur7 - 10thpctile (value of 25):
display a[21,21] // var of rur7
di a[30,30] // var of revXrur7
di a[21,30] //cov of rur7, revXrur7
di sqrt(a[21,21]+((25^2)*a[30,30])+(2*25*a[21,30])) // se=0.32325557
*LB CI - use non exponentiated coefficient
di exp((_b[rur7]+(_b[revXrur7]*25)) - 1.96*0.32325557) // =0.212
*UB CI - use non exponentiated coefficient
di exp((_b[rur7]+(_b[revXrur7]*25)) + 1.96*0.32325557) // =0.754

*revXrur7 - 50thpctile (value of 40):
di sqrt(a[21,21]+((40^2)*a[30,30])+(2*40*a[21,30])) // se=0.22639962
*LB CI - use non exponentiated coefficient
di exp((_b[rur7]+(_b[revXrur7]*40)) - 1.96*0.22639962) // =0.430
*UB CI - use non exponentiated coefficient
di exp((_b[rur7]+(_b[revXrur7]*40)) + 1.96*0.22639962) // =1.045

*revXrur7 - 90thpctile (value of 55):
di sqrt(a[21,21]+((55^2)*a[30,30])+(2*55*a[21,30])) // se=0.20619351
*LB CI - use non exponentiated coefficient
di exp((_b[rur7]+(_b[revXrur7]*55)) - 1.96*0.20619351) // =0.751
*UB CI - use non exponentiated coefficient
di exp((_b[rur7]+(_b[revXrur7]*55)) + 1.96*0.20619351) // =1.685





********** rur level 8 (8th most urban)
*revXrur8 - 10thpctile (value of 25):
di exp(_b[rur8]+(_b[revXrur8]*25)) //= 0.338
*revXrur8 - 50thpctile (value of 40):
di exp(_b[rur8]+(_b[revXrur8]*40)) //= 0.617
*revXrur8 - 90thpctile (value of 55):
di exp(_b[rur8]+(_b[revXrur8]*55)) //= 1.127

*revXrur8 - 10thpctile (value of 25):
display a[22,22] // var of rur8
di a[31,31] // var of revXrur8
di a[22,31] //cov of rur8, revXrur8
di sqrt(a[22,22]+((25^2)*a[31,31])+(2*25*a[22,31])) // se=0.35509152
*LB CI - use non exponentiated coefficient
di exp((_b[rur8]+(_b[revXrur8]*25)) - 1.96*0.35509152) // =0.168
*UB CI - use non exponentiated coefficient
di exp((_b[rur8]+(_b[revXrur8]*25)) + 1.96*0.35509152) // =0.677

*revXrur8 - 50thpctile (value of 40):
di sqrt(a[22,22]+((40^2)*a[31,31])+(2*40*a[22,31])) // se=0.25484913
*LB CI - use non exponentiated coefficient
di exp((_b[rur8]+(_b[revXrur8]*40)) - 1.96*0.25484913) // =0.374
*UB CI - use non exponentiated coefficient
di exp((_b[rur8]+(_b[revXrur8]*40)) + 1.96*0.25484913) // =1.017

*revXrur8 - 90thpctile (value of 55):
di sqrt(a[22,22]+((55^2)*a[31,31])+(2*55*a[22,31])) // se=0.22393535
*LB CI - use non exponentiated coefficient
di exp((_b[rur8]+(_b[revXrur8]*55)) - 1.96*0.22393535) // =0.727
*UB CI - use non exponentiated coefficient
di exp((_b[rur8]+(_b[revXrur8]*55)) + 1.96*0.22393535) // =1.748





********** rur level 9 (9th most urban)
*revXrur9 - 10thpctile (value of 25):
di exp(_b[rur9]+(_b[revXrur9]*25)) //= 0.222
*revXrur9 - 50thpctile (value of 40):
di exp(_b[rur9]+(_b[revXrur9]*40)) //= 0.471
*revXrur9 - 90thpctile (value of 55):
di exp(_b[rur9]+(_b[revXrur9]*55)) //= 0.998

*revXrur9 - 10thpctile (value of 25):
display a[23,23] // var of rur9
di a[32,32] // var of revXrur9
di a[23,32] //cov of rur9, revXrur9
di sqrt(a[23,23]+((25^2)*a[32,32])+(2*25*a[23,32])) // se=0.34942713
*LB CI - use non exponentiated coefficient
di exp((_b[rur9]+(_b[revXrur9]*25)) - 1.96*0.34942713) // =0.112
*UB CI - use non exponentiated coefficient
di exp((_b[rur9]+(_b[revXrur9]*25)) + 1.96*0.34942713) // =0.440

*revXrur9 - 50thpctile (value of 40):
di sqrt(a[23,23]+((40^2)*a[32,32])+(2*40*a[23,32])) // se=0.2503807
*LB CI - use non exponentiated coefficient
di exp((_b[rur9]+(_b[revXrur9]*40)) - 1.96*0.2503807) // =0.288
*UB CI - use non exponentiated coefficient
di exp((_b[rur9]+(_b[revXrur9]*40)) + 1.96*0.2503807) // =0.769

*revXrur9 - 90thpctile (value of 55):
di sqrt(a[23,23]+((55^2)*a[32,32])+(2*55*a[23,32])) // se=0.24562043
*LB CI - use non exponentiated coefficient
di exp((_b[rur9]+(_b[revXrur9]*55)) - 1.96*0.24562043) // =0.617
*UB CI - use non exponentiated coefficient
di exp((_b[rur9]+(_b[revXrur9]*55)) + 1.96*0.24562043) // =1.615





********** rur level 10 (10th most urban)
*revXrur10 - 10thpctile (value of 25):
di exp(_b[rur10]+(_b[revXrur10]*25)) //= 0.201
*revXrur10 - 50thpctile (value of 40):
di exp(_b[rur10]+(_b[revXrur10]*40)) //= 0.433
*revXrur10 - 90thpctile (value of 55):
di exp(_b[rur10]+(_b[revXrur10]*55)) //= 0.931

*revXrur10 - 10thpctile (value of 25):
display a[24,24] // var of rur10
di a[33,33] // var of revXrur10
di a[24,33] //cov of rur10, revXrur10
di sqrt(a[24,24]+((25^2)*a[33,33])+(2*25*a[24,33])) // se=0.39918225
*LB CI - use non exponentiated coefficient
di exp((_b[rur10]+(_b[revXrur10]*25)) - 1.96*0.39918225) // =0.092
*UB CI - use non exponentiated coefficient
di exp((_b[rur10]+(_b[revXrur10]*25)) + 1.96*0.39918225) // =0.440

*revXrur10 - 50thpctile (value of 40):
di sqrt(a[24,24]+((40^2)*a[33,33])+(2*40*a[24,33])) // se=0.26288312
*LB CI - use non exponentiated coefficient
di exp((_b[rur10]+(_b[revXrur10]*40)) - 1.96*0.26288312) // =0.258
*UB CI - use non exponentiated coefficient
di exp((_b[rur10]+(_b[revXrur10]*40)) + 1.96*0.26288312) // =0.724

*revXrur10 - 90thpctile (value of 55):
di sqrt(a[24,24]+((55^2)*a[33,33])+(2*55*a[24,33])) // se=0.34341479
*LB CI - use non exponentiated coefficient
di exp((_b[rur10]+(_b[revXrur10]*55)) - 1.96*0.34341479) // =0.475
*UB CI - use non exponentiated coefficient
di exp((_b[rur10]+(_b[revXrur10]*55)) + 1.96*0.34341479) // =1.825










***************************** NGOxRURAL INTERACTION ****************************

***************************** marginal effects *********************************

*Model 8 - ngoXpop_rural interaction
xtgee mob_tel i.inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 merch_ ///
rev_source fdi_pct_ expt_c ib2.igf_r c.ngo##c.rur_pop_pct ///
c.time##c.time if dev_sam==1, family(nbinomial 1.168) link(log) ///
corr(ar1) vce(robust) exposure(pop_tot) eform force

sum ngo if dev_s==1, d



*at 10th to 90th of rev source
margins, at(ngo=(167 499 1567) rur_pop_pct=(0(10)100) (means) _all) vsquish
mplotoffset, x(rur_pop_pct) offset(1.1) recastci(rspike) 



*qic value
gen ngoXrur = ngo*rur_pop_pct

qic mob_tel g2 g3 gdp_grow elect_r work_pop_chg fixed_tel_100 merch_ ///
fdi_pct_ expt_c rev_source igf1 igf2 rur_pop_pct ngo ngoXrur ///
time time2 if dev_sam==1, family(nbinomial 1.168) link(log) ///
corr(ar1) robust exposure(pop_tot) eform force i(id_cow) t(year) nodis





********************** Test using multiplicative effects ***********************

*Using procedure described by Hilbe (2011) pp.528-529

*costructing decileXngo interactions (1st decile is ref)
gen ngoXrur2 = rur2*ngo
gen ngoXrur3 = rur3*ngo
gen ngoXrur4 = rur4*ngo
gen ngoXrur5 = rur5*ngo
gen ngoXrur6 = rur6*ngo
gen ngoXrur7 = rur7*ngo
gen ngoXrur8 = rur8*ngo
gen ngoXrur9 = rur9*ngo
gen ngoXrur10 = rur10*ngo

*running model 8 with decile interactions
xtgee mob_tel i.inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 merch_ ///
fdi_pct_ expt_c rev_source ib2.igf_r ///
ngo rur2 rur3 rur4 rur5 rur6 rur7 rur8 rur9 rur10 ngoXrur2 ngoXrur3 ///
ngoXrur4 ngoXrur5 ngoXrur6 ngoXrur7 ngoXrur8 ngoXrur9 ngoXrur10 ///
c.time##c.time if dev_sam==1, family(nbinomial 1.168) link(log) ///
corr(ar1) vce(robust) exposure(pop_tot) eform force



*dropping svmat variables for next procedure
drop a1-a36



*list model vectors
ereturn list
matrix list e(b)

sum ngo if dev_s==1, d
***/// calculate multiplicative values for interactions 
***/// IRR=exp[b_dum+(b_int*ngo_value)]


***/// calculate CIs
***/// SE= sqrt(V(b_dum)+((ngo_val^2)*V(b_int))+(2*ngo_val*COV(b_dum,b_int)))
mat list e(V)
matrix a = e(V)
svmat a


**********////////// For figure 4 \\\\\\\\\\*********



********** rur level 2 (2nd most urban)
*ngoXrur2 - 10thpctile (value of 167):
di exp(_b[rur2]+(_b[ngoXrur2]*167)) //= 0.543
*ngoXrur2 - 50thpctile (value of 499):
di exp(_b[rur2]+(_b[ngoXrur2]*499)) //= 0.644
*ngoXrur2 - 90thpctile (value of 1567):
di exp(_b[rur2]+(_b[ngoXrur2]*1567)) //= 1.112

*ngoXrur2 - 10thpctile (value of 167):
display a[16,16] // var of rur2
di a[25,25] // var of ngoXrur2
di a[16,25] //cov of rur2, ngoXrur2
di sqrt(a[16,16]+((167^2)*a[25,25])+(2*167*a[16,25])) // se=0.23884137
*LB CI - use non exponentiated coefficient
di exp((_b[rur2]+(_b[ngoXrur2]*167)) - 1.96*0.23884137) // =0.340
*UB CI - use non exponentiated coefficient
di exp((_b[rur2]+(_b[ngoXrur2]*167)) + 1.96*0.23884137) // =0.867

*ngoXrur2 - 50thpctile (value of 499):
di sqrt(a[16,16]+((499^2)*a[25,25])+(2*499*a[16,25])) // se=0.19754185
*LB CI - use non exponentiated coefficient
di exp((_b[rur2]+(_b[ngoXrur2]*499)) - 1.96*0.19754185) // =0.437
*UB CI - use non exponentiated coefficient
di exp((_b[rur2]+(_b[ngoXrur2]*499)) + 1.96*0.19754185) // =0.948

*ngoXrur2 - 90thpctile (value of 1567):
di sqrt(a[16,16]+((1567^2)*a[25,25])+(2*1567*a[16,25])) // se=0.09361514
*LB CI - use non exponentiated coefficient
di exp((_b[rur2]+(_b[ngoXrur2]*1567)) - 1.96*0.09361514) // =0.925
*UB CI - use non exponentiated coefficient
di exp((_b[rur2]+(_b[ngoXrur2]*1567)) + 1.96*0.09361514) // =1.336





********** rur level 3 (3rd most urban)
*ngoXrur3 - 10thpctile (value of 167):
di exp(_b[rur3]+(_b[ngoXrur3]*167)) //= 0.494
*ngoXrur3 - 50thpctile (value of 499):
di exp(_b[rur3]+(_b[ngoXrur3]*499)) //= 0.621
*ngoXrur3 - 90thpctile (value of 1567):
di exp(_b[rur3]+(_b[ngoXrur3]*1567)) //= 1.291

*ngoXrur3 - 10thpctile (value of 167):
display a[17,17] // var of rur3
di a[26,26] // var of ngoXrur3
di a[17,26] //cov of rur3, ngoXrur3
di sqrt(a[17,17]+((167^2)*a[26,26])+(2*167*a[17,26])) // se=0.28076241
*LB CI - use non exponentiated coefficient
di exp((_b[rur3]+(_b[ngoXrur3]*167)) - 1.96*0.28076241) // =0.285
*UB CI - use non exponentiated coefficient
di exp((_b[rur3]+(_b[ngoXrur3]*167)) + 1.96*0.28076241) // =0.857

*ngoXrur3 - 50thpctile (value of 499):
di sqrt(a[17,17]+((499^2)*a[26,26])+(2*499*a[17,26])) // se=0.23388632
*LB CI - use non exponentiated coefficient
di exp((_b[rur3]+(_b[ngoXrur3]*499)) - 1.96*0.23388632) // =0.392
*UB CI - use non exponentiated coefficient
di exp((_b[rur3]+(_b[ngoXrur3]*499)) + 1.96*0.23388632) // =0.981

*ngoXrur3 - 90thpctile (value of 1567):
di sqrt(a[17,17]+((1567^2)*a[26,26])+(2*1567*a[17,26])) // se=0.120989
*LB CI - use non exponentiated coefficient
di exp((_b[rur3]+(_b[ngoXrur3]*1567)) - 1.96*0.120989) // =1.019
*UB CI - use non exponentiated coefficient
di exp((_b[rur3]+(_b[ngoXrur3]*1567)) + 1.96*0.120989) // =1.637





********** rur level 4 (4th most urban)
*ngoXrur4 - 10thpctile (value of 167):
di exp(_b[rur4]+(_b[ngoXrur4]*167)) //= 0.477
*ngoXrur4 - 50thpctile (value of 499):
di exp(_b[rur4]+(_b[ngoXrur4]*499)) //= 0.599
*ngoXrur4 - 90thpctile (value of 1567):
di exp(_b[rur4]+(_b[ngoXrur4]*1567)) //= 1.246

*ngoXrur4 - 10thpctile (value of 167):
display a[18,18] // var of rur4
di a[27,27] // var of ngoXrur4
di a[18,27] //cov of rur4, ngoXrur4
di sqrt(a[18,18]+((167^2)*a[27,27])+(2*167*a[18,27])) // se=0.30378841
*LB CI - use non exponentiated coefficient
di exp((_b[rur4]+(_b[ngoXrur4]*167)) - 1.96*0.30378841) // =0.263
*UB CI - use non exponentiated coefficient
di exp((_b[rur4]+(_b[ngoXrur4]*167)) + 1.96*0.30378841) // =0.865

*ngoXrur4 - 50thpctile (value of 499):
di sqrt(a[18,18]+((499^2)*a[27,27])+(2*499*a[18,27])) // se=0.25670629
*LB CI - use non exponentiated coefficient
di exp((_b[rur4]+(_b[ngoXrur4]*499)) - 1.96*0.25670629) // =0.362
*UB CI - use non exponentiated coefficient
di exp((_b[rur4]+(_b[ngoXrur4]*499)) + 1.96*0.25670629) // =0.991

*ngoXrur4 - 90thpctile (value of 1567):
di sqrt(a[18,18]+((1567^2)*a[27,27])+(2*1567*a[18,27])) // se=0.15615428
*LB CI - use non exponentiated coefficient
di exp((_b[rur4]+(_b[ngoXrur4]*1567)) - 1.96*0.15615428) // =0.918
*UB CI - use non exponentiated coefficient
di exp((_b[rur4]+(_b[ngoXrur4]*1567)) + 1.96*0.15615428) // =1.693





********** rur level 5 (5th most urban)
*ngoXrur5 - 10thpctile (value of 167):
di exp(_b[rur5]+(_b[ngoXrur5]*167)) //= 0.477
*ngoXrur5 - 50thpctile (value of 499):
di exp(_b[rur5]+(_b[ngoXrur5]*499)) //= 0.594
*ngoXrur5 - 90thpctile (value of 1567):
di exp(_b[rur5]+(_b[ngoXrur5]*1567)) //= 1.204

*ngoXrur5 - 10thpctile (value of 167):
display a[19,19] // var of rur5
di a[28,28] // var of ngoXrur5
di a[19,28] //cov of rur5, ngoXrur5
di sqrt(a[19,19]+((167^2)*a[28,28])+(2*167*a[19,28])) // se=0.3125031
*LB CI - use non exponentiated coefficient
di exp((_b[rur5]+(_b[ngoXrur5]*167)) - 1.96*0.3125031) // =0.259
*UB CI - use non exponentiated coefficient
di exp((_b[rur5]+(_b[ngoXrur5]*167)) + 1.96*0.3125031) // =0.881

*ngoXrur5 - 50thpctile (value of 499):
di sqrt(a[19,19]+((499^2)*a[28,28])+(2*499*a[19,28])) // se=0.26114726
*LB CI - use non exponentiated coefficient
di exp((_b[rur5]+(_b[ngoXrur5]*499)) - 1.96*0.26114726) // =0.356
*UB CI - use non exponentiated coefficient
di exp((_b[rur5]+(_b[ngoXrur5]*499)) + 1.96*0.26114726) // =0.992

*ngoXrur5 - 90thpctile (value of 1567):
di sqrt(a[19,19]+((1567^2)*a[28,28])+(2*1567*a[19,28])) // se=0.16635285
*LB CI - use non exponentiated coefficient
di exp((_b[rur5]+(_b[ngoXrur5]*1567)) - 1.96*0.16635285) // =0.869
*UB CI - use non exponentiated coefficient
di exp((_b[rur5]+(_b[ngoXrur5]*1567)) + 1.96*0.16635285) // =1.668





********** rur level 6 (6th most urban)
*ngoXrur6 - 10thpctile (value of 167):
di exp(_b[rur6]+(_b[ngoXrur6]*167)) //= 0.365
*ngoXrur6 - 50thpctile (value of 499):
di exp(_b[rur6]+(_b[ngoXrur6]*499)) //= 0.518
*ngoXrur6 - 90thpctile (value of 1567):
di exp(_b[rur6]+(_b[ngoXrur6]*1567)) //= 1.591

*ngoXrur6 - 10thpctile (value of 167):
display a[20,20] // var of rur6
di a[29,29] // var of ngoXrur6
di a[20,29] //cov of rur6, ngoXrur6
di sqrt(a[20,20]+((167^2)*a[29,29])+(2*167*a[20,29])) // se=0.32516047
*LB CI - use non exponentiated coefficient
di exp((_b[rur6]+(_b[ngoXrur6]*167)) - 1.96*0.32516047) // =0.193
*UB CI - use non exponentiated coefficient
di exp((_b[rur6]+(_b[ngoXrur6]*167)) + 1.96*0.32516047) // =0.691

*ngoXrur6 - 50thpctile (value of 499):
di sqrt(a[20,20]+((499^2)*a[29,29])+(2*499*a[20,29])) // se=0.27353672
*LB CI - use non exponentiated coefficient
di exp((_b[rur6]+(_b[ngoXrur6]*499)) - 1.96*0.27353672) // =0.303
*UB CI - use non exponentiated coefficient
di exp((_b[rur6]+(_b[ngoXrur6]*499)) + 1.96*0.27353672) // =0.885

*ngoXrur6 - 90thpctile (value of 1567):
di sqrt(a[20,20]+((1567^2)*a[29,29])+(2*1567*a[20,29])) // se=0.20479562
*LB CI - use non exponentiated coefficient
di exp((_b[rur6]+(_b[ngoXrur6]*1567)) - 1.96*0.20479562) // =1.065
*UB CI - use non exponentiated coefficient
di exp((_b[rur6]+(_b[ngoXrur6]*1567)) + 1.96*0.20479562) // =2.377





********** rur level 7 (7th most urban)
*ngoXrur7 - 10thpctile (value of 167):
di exp(_b[rur7]+(_b[ngoXrur7]*167)) //= 0.279
*ngoXrur7 - 50thpctile (value of 499):
di exp(_b[rur7]+(_b[ngoXrur7]*499)) //= 0.445
*ngoXrur7 - 90thpctile (value of 1567):
di exp(_b[rur7]+(_b[ngoXrur7]*1567)) //= 1.984

*ngoXrur7 - 10thpctile (value of 167):
display a[21,21] // var of rur7
di a[30,30] // var of ngoXrur7
di a[21,30] //cov of rur7, ngoXrur7
di sqrt(a[21,21]+((167^2)*a[30,30])+(2*167*a[21,30])) // se=0.32509827
*LB CI - use non exponentiated coefficient
di exp((_b[rur7]+(_b[ngoXrur7]*167)) - 1.96*0.32509827) // =0.148
*UB CI - use non exponentiated coefficient
di exp((_b[rur7]+(_b[ngoXrur7]*167)) + 1.96*0.32509827) // =0.528

*ngoXrur7 - 50thpctile (value of 499):
di sqrt(a[21,21]+((499^2)*a[30,30])+(2*499*a[21,30])) // se=0.27108505
*LB CI - use non exponentiated coefficient
di exp((_b[rur7]+(_b[ngoXrur7]*499)) - 1.96*0.27108505) // =0.261
*UB CI - use non exponentiated coefficient
di exp((_b[rur7]+(_b[ngoXrur7]*499)) + 1.96*0.27108505) // =0.756

*ngoXrur7 - 90thpctile (value of 1567):
di sqrt(a[21,21]+((1567^2)*a[30,30])+(2*1567*a[21,30])) // se=0.24013848
*LB CI - use non exponentiated coefficient
di exp((_b[rur7]+(_b[ngoXrur7]*1567)) - 1.96*0.24013848) // =1.239
*UB CI - use non exponentiated coefficient
di exp((_b[rur7]+(_b[ngoXrur7]*1567)) + 1.96*0.24013848) // =3.176





********** rur level 8 (8th most urban)
*ngoXrur8 - 10thpctile (value of 167):
di exp(_b[rur8]+(_b[ngoXrur8]*167)) //= 0.256
*ngoXrur8 - 50thpctile (value of 499):
di exp(_b[rur8]+(_b[ngoXrur8]*499)) //= 0.423
*ngoXrur8 - 90thpctile (value of 1567):
di exp(_b[rur8]+(_b[ngoXrur8]*1567)) //= 2.127

*ngoXrur8 - 10thpctile (value of 167):
display a[22,22] // var of rur8
di a[31,31] // var of ngoXrur8
di a[22,31] //cov of rur8, ngoXrur8
di sqrt(a[22,22]+((167^2)*a[31,31])+(2*167*a[22,31])) // se=0.3272913
*LB CI - use non exponentiated coefficient
di exp((_b[rur8]+(_b[ngoXrur8]*167)) - 1.96*0.3272913) // =0.135
*UB CI - use non exponentiated coefficient
di exp((_b[rur8]+(_b[ngoXrur8]*167)) + 1.96*0.3272913) // =0.487

*ngoXrur8 - 50thpctile (value of 499):
di sqrt(a[22,22]+((499^2)*a[31,31])+(2*499*a[22,31])) // se=0.28931715
*LB CI - use non exponentiated coefficient
di exp((_b[rur8]+(_b[ngoXrur8]*499)) - 1.96*0.28931715) // =0.240
*UB CI - use non exponentiated coefficient
di exp((_b[rur8]+(_b[ngoXrur8]*499)) + 1.96*0.28931715) // =0.747

*ngoXrur8 - 90thpctile (value of 1567):
di sqrt(a[22,22]+((1567^2)*a[31,31])+(2*1567*a[22,31])) // se=0.35016452
*LB CI - use non exponentiated coefficient
di exp((_b[rur8]+(_b[ngoXrur8]*1567)) - 1.96*0.35016452) // =1.071
*UB CI - use non exponentiated coefficient
di exp((_b[rur8]+(_b[ngoXrur8]*1567)) + 1.96*0.35016452) // =4.225





********** rur level 9 (9th most urban)
*ngoXrur9 - 10thpctile (value of 167):
di exp(_b[rur9]+(_b[ngoXrur9]*167)) //= 0.180
*ngoXrur9 - 50thpctile (value of 499):
di exp(_b[rur9]+(_b[ngoXrur9]*499)) //= 0.353
*ngoXrur9 - 90thpctile (value of 1567):
di exp(_b[rur9]+(_b[ngoXrur9]*1567)) //= 3.01

*ngoXrur9 - 10thpctile (value of 167):
display a[23,23] // var of rur9
di a[32,32] // var of ngoXrur9
di a[23,32] //cov of rur9, ngoXrur9
di sqrt(a[23,23]+((167^2)*a[32,32])+(2*167*a[23,32])) // se=0.33479178
*LB CI - use non exponentiated coefficient
di exp((_b[rur9]+(_b[ngoXrur9]*167)) - 1.96*0.33479178) // =0.093
*UB CI - use non exponentiated coefficient
di exp((_b[rur9]+(_b[ngoXrur9]*167)) + 1.96*0.33479178) // =0.347

*ngoXrur9 - 50thpctile (value of 499):
di sqrt(a[23,23]+((499^2)*a[32,32])+(2*499*a[23,32])) // se=0.29255792
*LB CI - use non exponentiated coefficient
di exp((_b[rur9]+(_b[ngoXrur9]*499)) - 1.96*0.29255792) // =0.199
*UB CI - use non exponentiated coefficient
di exp((_b[rur9]+(_b[ngoXrur9]*499)) + 1.96*0.29255792) // =0.627

*ngoXrur9 - 90thpctile (value of 1567):
di sqrt(a[23,23]+((1567^2)*a[32,32])+(2*1567*a[23,32])) // se=0.43106484
*LB CI - use non exponentiated coefficient
di exp((_b[rur9]+(_b[ngoXrur9]*1567)) - 1.96*0.43106484) // =1.332
*UB CI - use non exponentiated coefficient
di exp((_b[rur9]+(_b[ngoXrur9]*1567)) + 1.96*0.43106484) // =7.216





********** rur level 10 (10th most urban)
*ngoXrur10 - 10thpctile (value of 167):
di exp(_b[rur10]+(_b[ngoXrur10]*167)) //= 0.164
*ngoXrur10 - 50thpctile (value of 499):
di exp(_b[rur10]+(_b[ngoXrur10]*499)) //= 0.322
*ngoXrur10 - 90thpctile (value of 1567):
di exp(_b[rur10]+(_b[ngoXrur10]*1567)) //= 2.801

*ngoXrur10 - 10thpctile (value of 167):
display a[24,24] // var of rur10
di a[33,33] // var of ngoXrur10
di a[24,33] //cov of rur10, ngoXrur10
di sqrt(a[24,24]+((167^2)*a[33,33])+(2*167*a[24,33])) // se=0.36998161
*LB CI - use non exponentiated coefficient
di exp((_b[rur10]+(_b[ngoXrur10]*167)) - 1.96*0.36998161) // =0.079
*UB CI - use non exponentiated coefficient
di exp((_b[rur10]+(_b[ngoXrur10]*167)) + 1.96*0.36998161) // =0.339

*ngoXrur10 - 50thpctile (value of 499):
di sqrt(a[24,24]+((499^2)*a[33,33])+(2*499*a[24,33])) // se=0.29523455
*LB CI - use non exponentiated coefficient
di exp((_b[rur10]+(_b[ngoXrur10]*499)) - 1.96*0.29523455) // =0.180
*UB CI - use non exponentiated coefficient
di exp((_b[rur10]+(_b[ngoXrur10]*499)) + 1.96*0.29523455) // =0.573

*ngoXrur10 - 90thpctile (value of 1567):
di sqrt(a[24,24]+((1567^2)*a[33,33])+(2*1567*a[24,33])) // se=0.52119959
*LB CI - use non exponentiated coefficient
di exp((_b[rur10]+(_b[ngoXrur10]*1567)) - 1.96*0.52119959) // =1.008
*UB CI - use non exponentiated coefficient
di exp((_b[rur10]+(_b[ngoXrur10]*1567)) + 1.96*0.52119959) // =7.780









************************ robustness checks (appendix 2) ************************

*Mixed-effects negative binomial model - menbreg cannot specify time corr 
menbreg mob_tel i.inc_grp95 gdp_grow elect_r work_pop_chg ///
fixed_tel_100 merch_ ///
rev_source ngo fdi_pct_ expt_c ib2.igf_r rur_pop_pct ///
c.time##c.time if dev_sam==1, exposure(pop_tot) || id_:, ///
covariance(unst) vce(cluster id_) irr dispersion(mean)

testnl (_b[2.inc_grp95]=0) (_b[3.inc_grp95]=0)
testnl (_b[0.igf_recode]=0) (_b[1.igf_recode]=0)





*Unconditional NBFE
nbreg mob_tel i.id gdp_grow elect_r work_pop_chg ///
fixed_tel_100 merch_ ///
rev_source ngo rur_pop_pct fdi_pct_ expt_c ib2.igf_r ///
c.time##c.time if dev_sam==1, ///
vce(robust) exposure(pop_tot) irr

testnl (_b[0.igf_recode]=0) (_b[1.igf_recode]=0)





*subscriptions per person as DV (gaussian GEE)
gen mob_per_100 = (mob_tel/pop_tot)*100
gen pop_log = log(pop_tot)

*GEE model for gaussian DV
xtgee mob_per_100 i.inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 ///
merch_ rev_source fdi_pct_ expt_c ib2.igf_r ngo rur_pop_pct pop_log ///
time if inc_grp95!=. & inc_grp95!=4, family(gaussian) link(identity) ///
corr(ar1) vce(robust) force

test (_b[2.inc_grp95]=0) (_b[3.inc_grp95]=0)
test (_b[0.igf_recode]=0) (_b[1.igf_recode]=0)

*qic value
qic mob_per_100 g2 g3 gdp_grow elect_r work_pop_chg fixed_tel_100 merch_ ///
rev_source fdi_pct_ expt_c igf1 igf2 ngo rur_pop_pct pop_log ///
time time2 if inc_grp95!=. & inc_grp95!=4, family(gaussian) link(identity) ///
corr(ar1) robust eform force i(id_cow) t(year) nodis





*debt model to get alpha value for model below
nbreg mob_tel i.inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 merch_ ///
rev_source fdi_pct_ expt_c ib2.igf_r ngo rur_pop_pct debt_pct_ ///
c.time##c.time if inc_grp95!=. & inc_grp95!=4, ///
vce(robust) exposure(pop_tot) 

*GEE debt model
xtgee mob_tel i.inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 merch_ ///
rev_source fdi_pct_ expt_c ib2.igf_r ngo rur_pop_pct debt_pct_ ///
c.time##c.time if inc_grp95!=. & inc_grp95!=4, ///
family(nbinomial 1.248) link(log) ///
corr(ar1) vce(robust) exposure(pop_tot) eform force

testnl (_b[2.inc_grp95]=0) (_b[3.inc_grp95]=0)
testnl (_b[0.igf_recode]=0) (_b[1.igf_recode]=0)

*qic value
qic mob_tel g2 g3 gdp_grow elect_r work_pop_chg fixed_tel_100 merch_ ///
rev_source fdi_pct_ expt_c igf1 igf2 ngo rur_pop_pct debt_pct_ ///
time time2 if inc_grp95!=. & inc_grp95!=4, ///
family(nbinomial 1.248) link(log) ///
corr(ar1) robust exposure(pop_tot) eform force i(id_cow) t(year) nodis











************************ FDIxRURAL INTERACTION (APPENDIX 3)*********************

********************** Test using multiplicative effects ***********************

*Using procedure described by Hilbe (2011) pp.528-529

gen fdiXrur2 = rur2*fdi_pct_
gen fdiXrur3 = rur3*fdi_pct_
gen fdiXrur4 = rur4*fdi_pct_
gen fdiXrur5 = rur5*fdi_pct_
gen fdiXrur6 = rur6*fdi_pct_
gen fdiXrur7 = rur7*fdi_pct_
gen fdiXrur8 = rur8*fdi_pct_
gen fdiXrur9 = rur9*fdi_pct_
gen fdiXrur10 = rur10*fdi_pct_


xtgee mob_tel i.inc_grp95 gdp_grow elect_r work_pop_chg fixed_tel_100 merch_ ///
rev_source expt_c ngo ib2.igf_r ///
fdi_pct_ rur2 rur3 rur4 rur5 rur6 rur7 rur8 rur9 rur10 fdiXrur2 fdiXrur3 fdiXrur4 fdiXrur5 fdiXrur6 fdiXrur7 fdiXrur8 fdiXrur9 fdiXrur10 ///
c.time##c.time if dev_sam==1, family(nbinomial 1.168) link(log) ///
corr(ar1) vce(robust) exposure(pop_tot) eform force



*dropping svmat variables for next procedure
drop a1-a36



*list model vectors
ereturn list
matrix list e(b)

sum fdi_pct_ if dev_s==1, d
***/// calculate multiplicative values for interactions 
***/// IRR=exp[b_dum+(b_int*fdi_value)]

***/// calculate CIs
***/// SE= sqrt(V(b_dum)+((fdi_val^2)*V(b_int))+(2*fdi_val*COV(b_dum,b_int)))
mat list e(V)
matrix a = e(V)
svmat a



**********////////// For appendix 3 \\\\\\\\\\*********



********** rur level 2 (2nd most urban)
*fdiXrur2 - 10thpctile (value of 5):
di exp(_b[rur2]+(_b[fdiXrur2]*5)) //= 0.833
*fdiXrur2 - 50thpctile (value of 23):
di exp(_b[rur2]+(_b[fdiXrur2]*23)) //= 0.827
*fdiXrur2 - 91stpctile (value of 71):
di exp(_b[rur2]+(_b[fdiXrur2]*71)) //= 0.814

*fdiXrur2 - 10thpctile (value of 5):
display a[16,16] // var of rur2
di a[25,25] // var of fdiXrur2
di a[16,25] //cov of rur2, fdiXrur2
di sqrt(a[16,16]+((5^2)*a[25,25])+(2*5*a[16,25])) // se=0.14989374
*LB CI - use non exponentiated coefficient
di exp((_b[rur2]+(_b[fdiXrur2]*5)) - 1.96*0.14989374) // =0.621
*UB CI - use non exponentiated coefficient
di exp((_b[rur2]+(_b[fdiXrur2]*5)) + 1.96*0.14989374) // =1.117

*fdiXrur2 - 50thpctile (value of 23):
di sqrt(a[16,16]+((23^2)*a[25,25])+(2*23*a[16,25])) // se=0.14508565
*LB CI - use non exponentiated coefficient
di exp((_b[rur2]+(_b[fdiXrur2]*23)) - 1.96*0.14508565) // =0.623
*UB CI - use non exponentiated coefficient
di exp((_b[rur2]+(_b[fdiXrur2]*23)) + 1.96*0.14508565) // =1.10

*fdiXrur2 - 91stpctile (value of 71):
di sqrt(a[16,16]+((71^2)*a[25,25])+(2*71*a[16,25])) // se=0.15338483
*LB CI - use non exponentiated coefficient
di exp((_b[rur2]+(_b[fdiXrur2]*71)) - 1.96*0.15338483) // =0.603
*UB CI - use non exponentiated coefficient
di exp((_b[rur2]+(_b[fdiXrur2]*71)) + 1.96*0.15338483) // =1.10





********** rur level 3 (3nd most urban)
*fdiXrur3 - 10thpctile (value of 5):
di exp(_b[rur3]+(_b[fdiXrur3]*5)) //= 0.789
*fdiXrur3 - 50thpctile (value of 23):
di exp(_b[rur3]+(_b[fdiXrur3]*23)) //= 0.815
*fdiXrur3 – 91stpctile (value of 71):
di exp(_b[rur3]+(_b[fdiXrur3]*71)) //= 0.886

*fdiXrur3 - 10thpctile (value of 5):
display a[17,17] // var of rur3
di a[26,26] // var of fdiXrur3
di a[17,26] //cov of rur3, fdiXrur3
di sqrt(a[17,17]+((5^2)*a[26,26])+(2*5*a[17,26])) // se=0.19353664
*LB CI - use non exponentiated coefficient
di exp((_b[rur3]+(_b[fdiXrur3]*5)) - 1.96*0.19353664) // =0.54
*UB CI - use non exponentiated coefficient
di exp((_b[rur3]+(_b[fdiXrur3]*5)) + 1.96*0.19353664) // =1.154

*fdiXrur3 - 50thpctile (value of 23):
di sqrt(a[17,17]+((23^2)*a[26,26])+(2*23*a[17,26])) // se=0.18557603
*LB CI - use non exponentiated coefficient
di exp((_b[rur3]+(_b[fdiXrur3]*23)) - 1.96*0.18557603) // =0.566
*UB CI - use non exponentiated coefficient
di exp((_b[rur3]+(_b[fdiXrur3]*23)) + 1.96*0.18557603) // =1.172

*fdiXrur3 – 91stpctile (value of 71):
di sqrt(a[17,17]+((71^2)*a[26,26])+(2*71*a[17,26])) // se=0.18096015
*LB CI - use non exponentiated coefficient
di exp((_b[rur3]+(_b[fdiXrur3]*71)) - 1.96*0.18096015) // =0.621
*UB CI - use non exponentiated coefficient
di exp((_b[rur3]+(_b[fdiXrur3]*71)) + 1.96*0.18096015) // =1.263



********** rur level 4 (4th most urban)
*fdiXrur4 - 10thpctile (value of 5):
di exp(_b[rur4]+(_b[fdiXrur4]*5)) //= 0.701
*fdiXrur4 - 50thpctile (value of 23):
di exp(_b[rur4]+(_b[fdiXrur4]*23)) //= 0.752
*fdiXrur4 - 90thpctile (value of 71):
di exp(_b[rur4]+(_b[fdiXrur4]*71)) //= 0.906

*fdiXrur4 - 10thpctile (value of 5):
display a[18,18] // var of rur4
di a[27,27] // var of fdiXrur4
di a[18,27] //cov of rur4, fdiXrur4
di sqrt(a[18,18]+((5^2)*a[27,27])+(2*5*a[18,27])) // se=0.26747516
*LB CI - use non exponentiated coefficient
di exp((_b[rur4]+(_b[fdiXrur4]*5)) - 1.96*0.26747516) // =0.415
*UB CI - use non exponentiated coefficient
di exp((_b[rur4]+(_b[fdiXrur4]*5)) + 1.96*0.26747516) // =1.185

*fdiXrur4 - 50thpctile (value of 23):
di sqrt(a[18,18]+((23^2)*a[27,27])+(2*23*a[18,27])) // se=0.22798384
*LB CI - use non exponentiated coefficient
di exp((_b[rur4]+(_b[fdiXrur4]*23)) - 1.96*0.22798384) // =0.481
*UB CI - use non exponentiated coefficient
di exp((_b[rur4]+(_b[fdiXrur4]*23)) + 1.96*0.22798384) // =1.176

*fdiXrur4 – 91stpctile (value of 71):
di sqrt(a[18,18]+((71^2)*a[27,27])+(2*71*a[18,27])) // se=0.23086081
*LB CI - use non exponentiated coefficient
di exp((_b[rur4]+(_b[fdiXrur4]*71)) - 1.96*0.23086081) // =0.576
*UB CI - use non exponentiated coefficient
di exp((_b[rur4]+(_b[fdiXrur4]*71)) + 1.96*0.23086081) // =1.424





********** rur level 5 (5th most urban)
*fdiXrur5 - 10thpctile (value of 5):
di exp(_b[rur5]+(_b[fdiXrur5]*5)) //= 0.730
*fdiXrur5 - 50thpctile (value of 23):
di exp(_b[rur5]+(_b[fdiXrur5]*23)) //= 0.776
*fdiXrur5 - 90thpctile (value of 71):
di exp(_b[rur5]+(_b[fdiXrur5]*71)) //= 0.915

*fdiXrur5 - 10thpctile (value of 5):
display a[19,19] // var of rur5
di a[28,28] // var of fdiXrur5
di a[19,28] //cov of rur5, fdiXrur5
di sqrt(a[19,19]+((5^2)*a[28,28])+(2*5*a[19,28])) // se=0.22940942
*LB CI - use non exponentiated coefficient
di exp((_b[rur5]+(_b[fdiXrur5]*5)) - 1.96*0.22940942) // =0.465
*UB CI - use non exponentiated coefficient
di exp((_b[rur5]+(_b[fdiXrur5]*5)) + 1.96*0.22940942) // =1.144

*fdiXrur5 - 50thpctile (value of 23):
di sqrt(a[19,19]+((23^2)*a[28,28])+(2*23*a[19,28])) // se=0.21797321
*LB CI - use non exponentiated coefficient
di exp((_b[rur5]+(_b[fdiXrur5]*23)) - 1.96*0.21797321) // =0.506
*UB CI - use non exponentiated coefficient
di exp((_b[rur5]+(_b[fdiXrur5]*23)) + 1.96*0.21797321) // =1.190

*fdiXrur5 – 91stpctile (value of 71):
di sqrt(a[19,19]+((71^2)*a[28,28])+(2*71*a[19,28])) // se=0.20932899
*LB CI - use non exponentiated coefficient
di exp((_b[rur5]+(_b[fdiXrur5]*71)) - 1.96*0.20932899) // =0.607
*UB CI - use non exponentiated coefficient
di exp((_b[rur5]+(_b[fdiXrur5]*71)) + 1.96*0.20932899) // =1.380





********** rur level 6 (6th most urban)
*fdiXrur6 - 10thpctile (value of 5):
di exp(_b[rur6]+(_b[fdiXrur6]*5)) //= 0.704
*fdiXrur6 - 50thpctile (value of 23):
di exp(_b[rur6]+(_b[fdiXrur6]*23)) //= 0.707
*fdiXrur6 - 90thpctile (value of 71):
di exp(_b[rur6]+(_b[fdiXrur6]*71)) //= 0.716

*fdiXrur6 - 10thpctile (value of 5):
display a[20,20] // var of rur6
di a[29,29] // var of fdiXrur6
di a[20,29] //cov of rur6, fdiXrur6
di sqrt(a[20,20]+((5^2)*a[29,29])+(2*5*a[20,29])) // se=0.23355742
*LB CI - use non exponentiated coefficient
di exp((_b[rur6]+(_b[fdiXrur6]*5)) - 1.96*0.23355742) // =0.445
*UB CI - use non exponentiated coefficient
di exp((_b[rur6]+(_b[fdiXrur6]*5)) + 1.96*0.23355742) // =1.112

*fdiXrur6 - 50thpctile (value of 23):
di sqrt(a[20,20]+((23^2)*a[29,29])+(2*23*a[20,29])) // se=0.22963107
*LB CI - use non exponentiated coefficient
di exp((_b[rur6]+(_b[fdiXrur6]*23)) - 1.96*0.22963107) // =0.451
*UB CI - use non exponentiated coefficient
di exp((_b[rur6]+(_b[fdiXrur6]*23)) + 1.96*0.22963107) // =1.11

*fdiXrur6 – 91stpctile (value of 71):
di sqrt(a[20,20]+((71^2)*a[29,29])+(2*71*a[20,29])) // se=0.22126088
*LB CI - use non exponentiated coefficient
di exp((_b[rur6]+(_b[fdiXrur6]*71)) - 1.96*0.22126088) // =0.464
*UB CI - use non exponentiated coefficient
di exp((_b[rur6]+(_b[fdiXrur6]*71)) + 1.96*0.22126088) // =1.104







********** rur level 7 (7th most urban)
*fdiXrur7 - 10thpctile (value of 5):
di exp(_b[rur7]+(_b[fdiXrur7]*5)) //= 0.617
*fdiXrur7 - 50thpctile (value of 23):
di exp(_b[rur7]+(_b[fdiXrur7]*23)) //= 0.610
*fdiXrur7 - 90thpctile (value of 71):
di exp(_b[rur7]+(_b[fdiXrur7]*71)) //= 0.589

*fdiXrur7 - 10thpctile (value of 5):
display a[21,21] // var of rur7
di a[30,30] // var of fdiXrur7
di a[21,30] //cov of rur7, fdiXrur7
di sqrt(a[21,21]+((5^2)*a[30,30])+(2*5*a[21,30])) // se=0.23950737
*LB CI - use non exponentiated coefficient
di exp((_b[rur7]+(_b[fdiXrur7]*5)) - 1.96*0.23950737) // =0.386
*UB CI - use non exponentiated coefficient
di exp((_b[rur7]+(_b[fdiXrur7]*5)) + 1.96*0.23950737) // =0.987

*fdiXrur7 - 50thpctile (value of 23):
di sqrt(a[21,21]+((23^2)*a[30,30])+(2*23*a[21,30])) // se=0.23473617
*LB CI - use non exponentiated coefficient
di exp((_b[rur7]+(_b[fdiXrur7]*23)) - 1.96*0.23473617) // =0.385
*UB CI - use non exponentiated coefficient
di exp((_b[rur7]+(_b[fdiXrur7]*23)) + 1.96*0.23473617) // =0.966

*fdiXrur7 – 91stpctile (value of 71):
di sqrt(a[21,21]+((71^2)*a[30,30])+(2*71*a[21,30])) // se=0.23933847
*LB CI - use non exponentiated coefficient
di exp((_b[rur7]+(_b[fdiXrur7]*71)) - 1.96*0.23933847) // =0.369
*UB CI - use non exponentiated coefficient
di exp((_b[rur7]+(_b[fdiXrur7]*71)) + 1.96*0.23933847) // =0.942




********** rur level 8 (8th most urban)
*fdiXrur8 - 10thpctile (value of 5):
di exp(_b[rur8]+(_b[fdiXrur8]*5)) //= 0.554
*fdiXrur8 - 50thpctile (value of 23):
di exp(_b[rur8]+(_b[fdiXrur8]*23)) //= 0.564
*fdiXrur8 – 91stpctile (value of 71):
di exp(_b[rur8]+(_b[fdiXrur8]*71)) //= 0.592

*fdiXrur8 - 10thpctile (value of 5):
display a[22,22] // var of rur8
di a[31,31] // var of fdiXrur8
di a[22,31] //cov of rur8, fdiXrur8
di sqrt(a[22,22]+((5^2)*a[31,31])+(2*5*a[22,31])) // se=0.29249929
*LB CI - use non exponentiated coefficient
di exp((_b[rur8]+(_b[fdiXrur8]*5)) - 1.96*0.29249929) // =0.312
*UB CI - use non exponentiated coefficient
di exp((_b[rur8]+(_b[fdiXrur8]*5)) + 1.96*0.29249929) // =0.983

*fdiXrur8 - 50thpctile (value of 23):
di sqrt(a[22,22]+((23^2)*a[31,31])+(2*23*a[22,31])) // se=0.2709307
*LB CI - use non exponentiated coefficient
di exp((_b[rur8]+(_b[fdiXrur8]*23)) - 1.96*0.2709307) // =0.332
*UB CI - use non exponentiated coefficient
di exp((_b[rur8]+(_b[fdiXrur8]*23)) + 1.96*0.2709307) // =0.959

*fdiXrur8 – 91stpctile (value of 71):
di sqrt(a[22,22]+((71^2)*a[31,31])+(2*71*a[22,31])) // se=0.23410835
*LB CI - use non exponentiated coefficient
di exp((_b[rur8]+(_b[fdiXrur8]*71)) - 1.96*0.23410835) // =0.374
*UB CI - use non exponentiated coefficient
di exp((_b[rur8]+(_b[fdiXrur8]*71)) + 1.96*0.23410835) // =0.937





********** rur level 9 (9th most urban)
*fdiXrur9 - 10thpctile (value of 5):
di exp(_b[rur9]+(_b[fdiXrur9]*5)) //= 0.368
*fdiXrur9 - 50thpctile (value of 23):
di exp(_b[rur9]+(_b[fdiXrur9]*23)) //= 0.420
*fdiXrur9 – 91stpctile (value of 71):
di exp(_b[rur9]+(_b[fdiXrur9]*71)) //= 0.597

*fdiXrur9 - 10thpctile (value of 5):
display a[23,23] // var of rur9
di a[32,32] // var of fdiXrur9
di a[23,32] //cov of rur9, fdiXrur9
di sqrt(a[23,23]+((5^2)*a[32,32])+(2*5*a[23,32])) // se=0.26644659
*LB CI - use non exponentiated coefficient
di exp((_b[rur9]+(_b[fdiXrur9]*5)) - 1.96*0.26644659) // =0.218
*UB CI - use non exponentiated coefficient
di exp((_b[rur9]+(_b[fdiXrur9]*5)) + 1.96*0.26644659) // =0.621

*fdiXrur9 - 50thpctile (value of 23):
di sqrt(a[23,23]+((23^2)*a[32,32])+(2*23*a[23,32])) // se=0.25951696
*LB CI - use non exponentiated coefficient
di exp((_b[rur9]+(_b[fdiXrur9]*23)) - 1.96*0.25951696) // =0.253
*UB CI - use non exponentiated coefficient
di exp((_b[rur9]+(_b[fdiXrur9]*23)) + 1.96*0.25951696) // =0.699

*fdiXrur9 – 91stpctile (value of 71):
di sqrt(a[23,23]+((71^2)*a[32,32])+(2*71*a[23,32])) // se=0.26230868
*LB CI - use non exponentiated coefficient
di exp((_b[rur9]+(_b[fdiXrur9]*71)) - 1.96*0.26230868) // =0.357
*UB CI - use non exponentiated coefficient
di exp((_b[rur9]+(_b[fdiXrur9]*71)) + 1.96*0.26230868) // =1.00




********** rur level 10 (10th most urban)
*fdiXrur10 - 10thpctile (value of 5):
di exp(_b[rur10]+(_b[fdiXrur10]*5)) //= 0.313
*fdiXrur10 - 50thpctile (value of 23):
di exp(_b[rur10]+(_b[fdiXrur10]*23)) //= 0.367
*fdiXrur10 - 90thpctile (value of 71):
di exp(_b[rur10]+(_b[fdiXrur10]*71)) //= 0.562

*fdiXrur10 - 10thpctile (value of 5):
display a[24,24] // var of rur10
di a[33,33] // var of fdiXrur10
di a[24,33] //cov of rur10, fdiXrur10
di sqrt(a[24,24]+((5^2)*a[33,33])+(2*5*a[24,33])) // se=0.28144292
*LB CI - use non exponentiated coefficient
di exp((_b[rur10]+(_b[fdiXrur10]*5)) - 1.96*0.28144292) // =0.18
*UB CI - use non exponentiated coefficient
di exp((_b[rur10]+(_b[fdiXrur10]*5)) + 1.96*0.28144292) // =0.544

*fdiXrur10 - 50thpctile (value of 23):
di sqrt(a[24,24]+((23^2)*a[33,33])+(2*23*a[24,33])) // se=0.26932295
*LB CI - use non exponentiated coefficient
di exp((_b[rur10]+(_b[fdiXrur10]*23)) - 1.96*0.26932295) // =0.217
*UB CI - use non exponentiated coefficient
di exp((_b[rur10]+(_b[fdiXrur10]*23)) + 1.96*0.26932295) // =0.623

*fdiXrur10 – 91stpctile (value of 71):
di sqrt(a[24,24]+((71^2)*a[33,33])+(2*71*a[24,33])) // se=0.25832413
*LB CI - use non exponentiated coefficient
di exp((_b[rur10]+(_b[fdiXrur10]*71)) - 1.96*0.25832413) // =0.339
*UB CI - use non exponentiated coefficient
di exp((_b[rur10]+(_b[fdiXrur10]*71)) + 1.96*0.25832413) // =0.933











********************** Multiple Imputation (MICE) PROCEDURE ********************

*create debtXrurality interaction for below
gen rurXdebt = rur_pop_pct*debt_pct_

*open log to store MICE results (may take a while to run)
log using [file_name], replace


*examine missingness - debt measure is quite problematic
mdesc mob_tel inc_grp95 gdp_grow work_pop_chg fixed_tel_100 merch_ ///
rev_source fdi_pct_ expt_c igf_r ngo rur_pop_pct debt_pct_ ///
if inc_grp95!=. & inc_grp95!=4, any

mdesc mob_tel inc_grp95 gdp_grow work_pop_chg fixed_tel_100 merch_ ///
rev_source fdi_pct_ expt_c igf_r ngo rur_pop_pct debt_pct_ ///
if inc_grp95!=. & inc_grp95!=4, all

mdesc mob_tel inc_grp95 gdp_grow work_pop_chg fixed_tel_100 merch_ ///
rev_source fdi_pct_ expt_c igf_r ngo rur_pop_pct debt_pct_ ///
if inc_grp95!=. & inc_grp95!=4, none

mdesc mob_tel inc_grp95 gdp_grow work_pop_chg fixed_tel_100 merch_ ///
rev_source fdi_pct_ expt_c igf_r ngo rur_pop_pct debt_pct_ ///
if inc_grp95!=. & inc_grp95!=4



***/// declare variables for imputation
mi set flong
mi xtset id_cow year

mi register imputed rurXdebt gdp_grow work_pop_chg fixed_tel_100 merch_ ///
fdi_pct_ expt_c ngo rur_pop_pct debt_pct_ ext_debt_pct stat_cap ext_debt_dol ///
polity2 hdi trade_pct rev_source



***/// Perform MICE (using predictive mean matching-to account for non-continuous variables and a neighbor pool of 6 for each draw)
mi impute chained (pmm, knn(6)) gdp_grow rurXdebt work_pop_chg ///
fixed_tel_100 merch_ fdi_pct_ expt_c ngo rur_pop_pct debt_pct_ ext_debt_pct ///
stat_cap ext_debt_dol polity2 hdi trade_pct rev_source = ///
mob_tel ib2.igf_r elect_r cs_repress urb_rur_excl ///
fh_civlib fh_polright igo gdp_pc pop_tot ///
i.inc_grp i.inc_grp95 time time2 i.id i.year ///
if inc_grp95!=. & inc_grp95!=4, ///
rseed(2182082) burnin(30) add(60) force dots

mi describe



*****///// MICE MODELS \\\\\*****

*model 6 from table 2
mi estimate, errorok eform saving(estfile1, replace) esample(esample531) ufmitest: xtgee mob_tel i.inc_grp95 gdp_grow elect_r work_pop_chg ///
fixed_tel_100 merch_ rev_source fdi_pct_ expt_c ib2.igf_r ngo ///
rur_pop_pct time time2 if inc_grp95!=. & inc_grp95!=4, ///
family(nbinomial) link(log) ///
corr(ar1) vce(robust) exposure(pop_tot) force

*Von Hippels command, indicates more than enough imputations
how_many_imputations

mi test 2.inc_grp95 3.inc_grp95, ufmitest
mi test 0.igf_recode 1.igf_recode, ufmitest

mi estimate (irr_incg2: exp(_b[2.inc_grp95])-1) using estfile1, ///
ufmitest
mi estimate (irr_incg3: exp(_b[3.inc_grp95])-1) using estfile1, ///
ufmitest

*check diagnostics
mi estimate, vartable dftable





*with Debt measure
mi estimate, errorok eform saving(estfile2, replace) esample(esample231) ufmitest: xtgee mob_tel i.inc_grp95 gdp_grow elect_r work_pop_chg ///
fixed_tel_100 merch_ rev_source fdi_pct_ expt_c ib2.igf_r ///
ngo rur_pop_pct debt_pct_ time time2 if inc_grp95!=. & inc_grp95!=4, ///
family(nbinomial) link(log) ///
corr(ar1) vce(robust) exposure(pop_tot) force

*Von Hippels command, enough imputations
how_many_imputations

mi test 2.inc_grp95 3.inc_grp95, ufmitest
mi test 0.igf_recode 1.igf_recode, ufmitest

mi estimate (irr_incg2: exp(_b[2.inc_grp95])-1) using estfile2, ///
ufmitest
mi estimate (irr_incg3: exp(_b[3.inc_grp95])-1) using estfile2, ///
ufmitest

*check diagnostics
mi estimate, vartable dftable





*with debt and rurality interaction (JAV approach)
mi estimate, errorok eform saving(estfile2, replace) esample(esample232) ufmitest: xtgee mob_tel i.inc_grp95 gdp_grow elect_r work_pop_chg ///
fixed_tel_100 merch_ rev_source fdi_pct_ expt_c ib2.igf_r ///
ngo rur_pop_pct debt_pct_ rurXdebt time time2 if inc_grp95!=. & inc_grp95!=4, ///
family(nbinomial) link(log) ///
corr(ar1) vce(robust) exposure(pop_tot) force

*Von Hippels command, enough imputations
how_many_imputations

mi test 2.inc_grp95 3.inc_grp95, ufmitest
mi test 0.igf_recode 1.igf_recode, ufmitest

mi estimate (irr_incg2: exp(_b[2.inc_grp95])-1) using estfile2, ///
ufmitest
mi estimate (irr_incg3: exp(_b[3.inc_grp95])-1) using estfile2, ///
ufmitest

log close


